Efficient measurement of linear susceptibilities in molecular 
simulations: Application to aging supercooled liquids 
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We propose a new method to measure time-dependent linear susceptibilities in molecular sim- 
ulations, which does not require the use of nonequilibrium simulations, subtraction techniques, or 
fluctuation-dissipation theorems. The main idea is an exact reformulation of linearly perturbed 
quantities in terms of observables accessible in unperturbed trajectories. We have applied these 
ideas to two supercooled liquids in their nonequilibrium aging regime. We show that previous work 
had underestimated deviations from fluctuation-dissipation relations in the case of a Lennard- Jones 
system, while our results for silica are in qualitative disagreement with earlier results. 
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Correlation and response functions play a major role 
in condensed matter physics as they directly probe static 
and dynamic properties at a microscopic level [lj. At 
thermal equilibrium, linear response theory permits the 
derivation of fluctuation-dissipation relations between 
conjugated susceptibilities and correlations, so that both 
types of measurements become equivalent [2[ . Depending 
on the technique used, experiments or simulations access 
one or the other quantity. For liquids, neutron scattering 
experiments will for instance be sensitive to spontaneous 
fluctuations of the density, while dielectric spectroscopy 
detects the response induced by an electric field 0]. Nu- 
merical simulations mainly focus on spontaneous fluctu- 
ations andprobe microscopic dynamics via correlation 
functions [3(. However, there exist cases where the nu- 
merical measurement of response functions becomes nec- 
essary, for instance when correlation functions become 
too noisy to be detected [H], or in nonequilibrium sit- 
uations, where correlation and response functions con- 
tain distinct physical information because fluctuation- 
dissipation theorems (FDT) do not hold [5[ . Quantifying 
FDT "violations" from the simultaneous measurement of 
correlation and response functions is an active field of re- 
search Q . In this work we propose an efficient method to 
access linear response functions in numerical simulations 
of molecular systems. As a physically relevant situation 
we apply this novel technique to study response functions 
of glass-forming liquids undergoing physical aging after 
a sudden quench to low temperature. 

Direct measurements of linear susceptibilities usually 
proceed as follows. Consider of system of N particles de- 
scribed by coordinates, r= {fj, i = 1, • • • , N}, momenta, 
p = {pi,i = 1, • • • , N}, masses m,, and a Hamiltonian 
H(r,p) containing a kinetic part, fC(p) = ^ZiVi /(^ m i)^ 
and a potential part, V(r). We first consider Newtonian 
dynamics, as used in Molecular Dynamics (MD): 

k = dH/dpi, # = -dH/dqi. (1) 
Physical observables, A(t) = A\p(t), f(t)], can be mea- 



sured at any time in a simulation, and correlation func- 
tions, C(t,t') = (A(t)B(t'))o, are obtained by averaging 
over repeated measurements. The subscript "0" indi- 
cates that averages are performed over unperturbed tra- 
jectories, and we suppose that (A(t))o = 0. In systems 
which are time-translationally invariant, two-time quan- 
tities only depend on t — t' but we retain the (i, t') nota- 
tion as we shall also study non-stationary systems. 

To measure a response function, an external field of 
constant amplitude h, conjugated to B(t), is introduced 
at time t' , such that the Hamiltonian contains the addi- 
tional term SH = —hB for t > t'. A linear susceptibility 
can then be defined: 



x(M') 



dh(t" 



(2) 



h->0 



Step responses are considered for simplicity but the dis- 
cussion holds more generally. The average in @ is with 
the field switched on, the zero-field limit comes from re- 
peated measurements with fields of decreasing amplitude. 
In practice, a compromise is sought between large fields 
introducing unwanted non-linear effects, and small fields 
resulting in poor signals. Such a non-equilibrium tech- 
nique suffers from a serious drawback. Averages in ^ are 
taken over perturbed trajectories, so that susceptibilities 
can only be measured one at a time, contrary to corre- 
lation functions which can be simultaneously measured 
and time averaged in a single unperturbed trajectory. 

An alternative would be to perform the derivative in 
Eq. @ before taking the average. This is precisely how 
the FDT is derived Averages are first expressed in 
terms of the distribution function. Its thermal equilib- 
rium (Gibbs-Boltzmann) form at temperature T is then 
assumed, and the derivative is computed analytically : 



x{t,t') = ^[c{t,t)-c{t,t')\ 



(3) 



where we have set Boltzmann's constant to unity. An 
important and well-known feature of the FDT in Eq. (J3]) 
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is that the right hand side is evaluated using unperturbed 
trajectories, the temperature prefactor reminding us that 
thermal equilibrium is assumed, implying that Eq. © 
cannot be used to measure far from equilibrium. 

The idea introduced in this paper is to perform the 
derivative before doing the average without assuming 
thermal equilibrium. Similar ideas were recently dis- 
cussed for discrete spins 0- In MD simulations, the 
subtraction technique is a finite-field approximation 
of this idea: Two simulations are run in parallel starting 
from the same configuration at time t' , one with h — 0, 
the other with a small field, h. The susceptibility reads: 
X(M') ~ {(A(t)) - (A(t)) ) /h. Non-equilibrium tech- 
niques are in fact unnecessary Q, since the h — > limit 
can be taken directly from (TT]) using perturbation the- 
ory Q to devise an unperturbed technique. The quanti- 
ties Xi = 9fi I dh and tfi = dpi / dh evolve as [§] : 

_$i dB(r,p) ^ _ dB(r,p) A d 2 V{r) _ 
Xl mi dp, ' % dr t 4-f df.drj Xr 

(4) 

The susceptibility x(^ ; t') can now be evaluated from un- 
perturbed trajectories: 

^)-(f(5ffi-* + ^.*)) o . < 5 » 

To illustrate the result in Eq. ([5]) we have performed 
MD simulations of a 80:20 binary Lennard- Jones (LJ) 
system composed of N = 10 3 particles at density p = 1.2. 
Particles interact with a LJ potential with parameters 
that can be found in [l(| , chosen to avoid crystallization 
at low temperature, and to study the properties of glass- 
forming liquids. Technical details of our simulations are 
as in the original paper [n| ■ When the temperature gets 
lower than T w 1 (we use LJ units the dynam- 

ics dramatically slows down, and the system cannot be 
equilibrated in computer simulations below T sa 0.43. 

We perform equilibrium simulations where we simul- 
taneously solve ((T|) and (H} to evaluate x{t,t') from (J5J), 
and the correlation C(t,t'). We focus on the follow- 
ing observables: A(t) = N^ 1 J^V ; ej cxp[ik ■ fj(t)] and 

B{t) = 2 J2j e j cos[fc • fj(t)], where ej = ±1 is a bimodal 
random variable of mean 0, such that C{t,t') corre- 
sponds to the self- intermediate scattering function 0. 
The numerical burden is a mere factor two since one in- 
tegrates 12iV instead of 6N equations of motion. For 
T = 1.0, dynamics is fast and X can b e evaluated 
in a few runs, as can be checked using the FDT. For 
T = 0.75, where the relaxation time is w 50 (see inset 
of Fig. [l}, the fundamental limitation of the technique 
appears. In Fig. [T] we represent T\{t,t') evaluated from 
10 3 independent runs using ((5]), as a function of C(t,t'). 
FDT predicts the linear relation shown as a full line. For 
t - t' < 5, x(t, t') follows the FDT. For larger t - t' , the 




C(t,t') 

FIG. 1: Simultaneous measurement of susceptibility \(t,t') 
and correlation C(t,t') in 10 3 independent unperturbed tra- 
jectories at T — 0.75 in the LJ system using Eq. ((SJ) for MD 
and Eq. ((6)1 for MC. For MD the noise diverges exponentially 
and x(M') cannot be evaluated for t — t' > 10, as indicated 
in the inset showing C(t,t') measured in MD. In MC simu- 
lations x(*>^') perfectly follows the FDT prediction indicated 
by a full line over the whole time range. 

noise in the susceptibility diverges exponentially and no 
reliable measurement can be performed, as in subtrac- 
tion techniques. Because the system is chaotic, nearby 
trajectories diverge exponentially quickly: While linear 
response fails at the level of trajectories [ll[, it holds 
at the probabilistic level Q, as suggested by the FDT 
derivation outlined above. 

The above exercise suggests that in Monte-Carlo (MC) 
simulations, where phase space is sampled probabilis- 
tically rather than deterministically, response functions 
could be efficiently evaluated. In a standard MC sim- 
ulation a configuration, Ct, is reached at time t. 
A trial configuration, C' t , is accessed with acceptance 
rate Ac t — >c' i generally defined from the energy change 
between the two configurations, e.g. the Metropolis 
rule 0] used in the following. The transition probabil- 
ity from Ct to C t+ i reads: W Ct ->c t+ i = ^C t+1 ,C' t A Ct ^c' t + 
<$c t+1 .c t (l — Ac t ^c^)- Averages now mean sampling 
a large number, N, of trajectories, (A(t)B(t'))o = 

JV _1 E£ii4k(*)£fc(tW*' -» *), where Mt) is the 
value of A at time t in trajectory k, and Pk(t — » t') is the 
probability of trajectory k between times t' and t starting 
from C t >, P k {t' ->■ t) = ]lt'~=t' Wc t fe „^e*„ ^ wherc C t" is 
the configuration visited at time t" in trajectory k. The 
susceptibility reads x(M') = ^[A/ --1 Y^k A k {t)P k {t' -> 
t)], and can be reformulated as an unperturbed average, 

X (t,t') = (A(t)R(t' ^t)) , (6) 

where R(t' ->■ t) = £ 4 „ d h M W C fe, ^ c \„ 1 ) • In Fig. CD 
we report the simultaneous measurement of x(^i*')> es ~ 
timated via ((6]), and of C(t,t') using 10 3 independent 
MC runs of the binary Lennard-Jones mixture described 
above for T = 0.75. (The details of the numerics ap- 
peared recently (l2l|.) The measurement now easily ex- 
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tends over the whole range of timescale over which C(t, t') 
changes, and FDT is perfectly obeyed. Although MC 
trajectories are chaotic, no exponential divergence of the 
noise is observed, at variance with the MD case. What 
Eq. © in fact does is to use a single unperturbed tra- 
jectories to evaluate the value the observable A(t) would 
have taken if an infinitesimal field had been applied. Ad- 
ditionally, the evaluation of Eq. (J6j) is computationally 
free since it only requires updating one additional observ- 
able, R(t' — > t), during the production of unperturbed 
trajectories. Finally, several susceptibilities and corre- 
lations may now be computed during the same simula- 
tion, and time averaging is easily implemented. The main 
limitation of the method is again statistics: x(t, t 1 ) now 
takes the form of a multi-time correlator, and its mea- 
surement becomes statistically costly as t — t' gets too 
large. We find an algebraic growth of the noise, as in spin 
systems 0, which is nevertheless a drastic improvement 
over exponential growth. A second drawback is the need 
to replace Newtonian by Monte-Carlo dynamics since the 
resulting dynamics are not necessarily equivalent. Quan- 
titative agreement between MC and MD dynamics was 
recently reported for the L J system described above [l2| • 
We now apply Eq. {BJ to measure x{t, t') after a sudden 
quench to very low temperature. Physical properties of 
the system now depend on the time t' spent since the 
quench, the system "ages" Energy slowly decreases 
with time, while dynamics gets slower [li|. The FDT in 
Eq. ([3]) no more applies, and the following generalization 



was suggested for glassy materials [14 1 



X{t,t') d 
T dt 



7 C(t,t'), 



(7) 



where X(t, t') is the fluctuation-dissipation ratio (FDR), 
X(t,t') — 1 at equilibrium. Deviations of the FDR from 
unity serve to quantify the distance from equilibrium [lI ] . 

Earlier attempts to measure X(t, t') in molecular 
glasses 15, 1(| used the following protocol: quench the 
system at t' = 0; apply a small field and measure t') 
for times t > t'; build a parametric "FD plot" of x{t,t') 
vs C(t, £'). Crucially, this amounts to replacing d t > by d t 
in ([7]), a procedure which is correct if X(t,t') is not an 
explicit function of t and t' [13] ■ Unbiased FDR measure- 
ments require instead the evaluation of x(£)0 & t fixed 
time t for different t' , so that the FDR can be graphi- 
cally deduced from the slope, —X(t,t')/T, of FD plots. 
This is numerically too costly if non-equilibrium tech- 
niques are used. The difficulty is easily overcome with 
Eq. ©, and we shall therefore report the first unbiased 
FDR measurements in aging molecular liquids. 

In Fig. [2] we use both time parametrizations to build 
FD plots in two glass-formers: the LJ system described 
above, and the BKS model for silica 1181 . The LJ results 




are qualitatively consistent with earlier reports 15j . The 
plots consist of two distinct pieces, FDT being satisfied 
for small t — t', "violated" for large t — t'. Strikingly, FD 



FIG. 2: Simultaneous measurement of x(*iO an d C(t,t') 
in aging LJ (T = 0.4, k = 6.7) and silica (T = 2500 K, 
k — 2.7 A -1 ). Fitting the nonequilibrium part of the FD 
plots (dashed line) for fixed-i parametrizations directly yields 
the FDRs x = 0.29 (LJ) and x = 0.49 (BKS). Incorrectly 
extracting x from fixed-i' data would yield 0.36 (LJ) and 0.63 
(BKS), seriously underestimating FDT deviations. 



plots are well-described by two straight lines, leading to 
a sensible definition of a constant FDR, x, at large t — t'. 
However, it is obvious in Fig. [2] that (incorrectly) esti- 
mating x from fixed-t' measurements yields values that 
seriously differ from unbiased estimates from fixed-i data, 
an error made in all previous FDR measurements [l5| . 
Both estimates only become equivalent if a non-trivial 
limiting FD plot is found at large time 14 1. 

For silica, we find similar FD plots, and similar quanti- 
tative discrepancies between both time parametrizations. 
The disagreement with earlier results is more pronounced 
for silica since FDR larger than unity were reported [n| . 
We have repeated our measurements at several tempera- 
tures between 500 and 2500 K, several wavevectors from 
0.3 to 13 A -1 , both for Si and O atoms. We consis- 
tently find FD plots as in Fig. [2] with X(t,t') < 1. We 
have numerically checked that this discrepancy cannot 
be explained by non-linear effects potentially present in 
the data of Ref. [lji|. Using non-equilibrium techniques 
with large fields we find that non-linear effects yield even 
smaller apparent FDR values. 

We have used the flexibility offered by Eq. © to char- 
acterize further the properties of FDRs in both aging 
liquids in Fig. [3J The top panel presents evidence that 
different observables share the same FDR value, obtained 
by changing the wavevector used to evaluate dynamic 
functions. Similar results were obtained for silica. The 
middle panel shows that Si and O atoms in silica display 
similar FD plots, with equal FDR values. Again, we find 
similar results for the two types of particles in the LJ 
mixture. These results suggest that it is sensible to de- 
fine, for fixed t, a unique FDR value x(t) characterizing 
the non-equilibrium part of FD plots. These findings are 
therefore compatible with the physical idea [l| that slow 
rearrangements in aging supercooled liquids behave as if 
they were thcrmalized at an "effective temperature" de- 
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FIG. 3: Top: FD plots for fixed T and t in the LJ system 
and different wavevectors displaying the same nonequilibrium 
value of the FDR. Middle: FD plots for Si and O (horizontally 
shifted by 0.2) in BKS for fixed T,k = 2.7k" 1 , and various t. 
For t = 4.10 4 , the FDR x = 0.51 fits both sets of data. Bot- 
tom: Temperature dependence of the FDR at a single large 
time, x(t = 10 4 ), for LJ and BKS systems. The temperature 
is normalized by the mode-coupling temperature T c . A linear 
behaviour (dashed line) is observed at low T. 



fined by T eS (t) = T/x{t) with T cff (i) > T in the 
two investigated systems. Our data indicate that T c g(t) 
decreases very slowly with t. Finally, the bottom panel 
shows the temperature dependence of the FDR measured 
at a single large time, x(t — 10 4 ). To compare both liq- 
uids we have to normalize the temperature by some tem- 
perature scale. We choose the "mode-coupling" temper- 
ature [T c = 0.435 (LJ) and T c = 3300 K (BKS)] because 
equilibration is numerically difficult below T c and aging 
effects can be detected. Remarkably, we find that FDRs 
in the two liquids display a very similar temperature de- 
pendence, x fss 0A7T/T C at small T. This confirms that 
both fragile (LJ) and strong (BKS silica) glass-formers 



studied in this work display similar aging properties. 

We have introduced a new technique to efficiently 
measure linear susceptibilities in molecular simulations 
which only uses unperturbed trajectories to evaluate re- 
sponse functions and outperforms subtraction techniques 
in Monte-Carlo simulations. Applied to aging super- 
cooled liquids, the technique allowed us to report the first 
unbiased numerical estimates of FDRs in aging molecular 
liquids, and to extend its determination to a wide range 
of times, temperatures, and observables. We showed that 
previous analysis quantitatively underestimated FDT vi- 
olations in LJ systems, while our results for silica are in 
qualitative disagreements with earlier results. 

I thank J.-L. Barrat who suggested to reconsider the 
aging regime of BKS silica and followed this work, and 
R.L. Jack and W. Kob for useful discussions and remarks 
on the manuscript. 
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